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ABSTRACT 

We present Particle-Based Lensing (PBL), a new technique for gravitational lensing mass recon- 
structions of galaxy clusters. Traditionally, most methods have employed either a finite inversion or 
gridding to turn observational lensed galaxy ellipticities into an estimate of the surface mass density 
of a galaxy cluster. We approach the problem from a different perspective, motivated by the success of 
multi-scale analysis in smoothed particle hydrodynamics. In PBL, we treat each of the lensed galaxies 
as a particle and then reconstruct the potential by smoothing over a local kernel with variable smooth- 
ing scale. In this way, we can tune a reconstruction to produce constant signal-to-noise throughout, 
and maximally exploit regions of high information density. 

PBL is designed to include all lensing observables, including multiple image positions and fluxes from 
strong lensing, as well as weak lensing signals including shear and flexion. In this paper, however, we 
describe a shear-only reconstruction, and apply the method to several test cases, including simulated 
lensing clusters, as well as the well-studied "Bullet Cluster" (1E0657-56). In the former cases, we 
show that PBL is better able to identify cusps and substructures than are grid-based reconstructions, 
and in the latter case, we show that PBL is able to identify substructure in the Bullet Cluster without 
even exploiting strong lensing measurements. We also make our codes publicly available. 
Subject headings: gravitational lensing, galaxiesxlusters 
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1. INTRODUCTION 

Cluster s of galaxies are excellent cosmological labo- 
rator ies (jAllen et all 120071 : iKingl l2007t iMannucci et all 
l2007f ). For example, the mass function of clus- 
ters is a sensitiv e probe of cosmo logical parameters 
like Q. m and <r s (Ri nes et al.1 |2007) and its observed 
evolution is an important test of theories of struc- 
ture formation dGunn fe GottJ 119721: iGiqcoli et alj 120071: 
IHorellou fe Bergdl2005l: ICoorav fc Shetbll2002D . The ge- 
ometrical shape of cluster Dark Matter halos provide 
valuable informat i on on intra-cluster gas distribution 
(|Flores et al.ll2005i r2007). While simulations predict cen- 
tral density distribution of matter in clusters to follow an 
NFW profile, it is debatable whethe r observations su, 
gests that cluste r s have a central core (JSand et al.l ( 200 
IVoigt fc Fabian] (|200l h 

ACDM structure formation theories also predict that 
massive dark matter halos assemble from the hierarchical 
merging of lower mass subh a los. As noted by s everal 
authors ([Moore et alj (jl999D : iKlypin et ail (|1999D h the 
number of subhalos that survive in N-body simulations is 
much greater than the number of dwarf galaxies observed 
in the Milky Way and the Andromeda galaxy. On cluster 
scales, such discrepancies are not observed. Thus the 
subhalo mass function in clusters is an important probe 
of the CDM theory in this mass scale. 

High-resolution, accurate measurements of cluster 
mass maps are thus highly desirable. Gravitational 
lensing is a powerful tool to probe the projected 
mass map of the clusters independent of the inter- 
nal dynamics, and has already been wi dely applied to 
mapping mass distribution in clusters ( Wittma n et al 



2004 iBroadhurst et all 120051: [Leonard et all 120071 : 
Okura et al.l 120071: iHevmans et all I2008T). Some re- 



120011 [Hockstr aet alj200 H iGrav et al.ll2002l : lTavlor etal 



search ers ( Nataraian fc Springe] l2004t Nataraian et al.1 
l2007bf l have used the individual galaxy-galaxy lensing 
signal to estimate individual galaxy masses and thus 
produce a parametric mass reconstruction of the clus- 
ter. Others have used the weak signal to characterize 
the overall potential from the cluster without recourse 
to parametric models dWilsori et al.lll996tlHoekstra et al.1 
Il998t iNataraian fc Refregierj|2000UHoekstra et alj2 004) . 

Given the importance of accurately measuring the 
mass, shape, and substructure of individual clusters, and 
given the enormous expense of long time-exposure ob- 
servations of clusters, it is extremely important to max- 
imize the signal-to-noise from a particular dataset, and 
to produce high-resolution maps of substructure within 
individual clusters. Current mass reconstruction tech- 
niques are ill-equipped to handle multi-scale datasets or 
clusters with significant dumpiness or cuspiness, or are 
jury-rigged to do so. In this paper, we propose Particle 
Based Lensing (PBL; pronounced "pebble") as an alter- 
native approach to cluster reconstruction. 

Our outline is as follows. In §[^]we give a brief review of 
the essential lensing formalism, and lay out our notation 
for the rest of the work. In§[3]we describe current (grid- 
based) techniques for reconstructing galaxy clusters, and 
identify some strengths and complications. In § [4] we 
propose Particle Based Lensing. We then apply this new 
method to simple simulated clusters of single and double 
peak softened isothermal spheres and the "bullet cluster" 
(1E0657-56) in § [5] We conclude in § [6] with a discus- 
sion of future prospects, including how additional strong- 
lensing and flexion information can be incorporated into 
PBL. 



Electronic address: sd365@drexel.edu 
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2. BACKGROUND 

Before delving into technical details of our method 
we would like to introduce the basic lensing nota- 
tion to be used throughout the paper. Following 
iBartelmann fc SchneideiT(|200lD . we consider a surface 
mass density £(0), where 9 is the angular position in the 
lens plane. Convergence or dimensionless surface mass 
density is defined as 



k{9) = 



where 



£(<?) 



(1) 



(2) 



AttG D d D ds 

D s , Dei, and D^s are the angular diameter distances be- 
tween the observer and the source, the observer and the 
lens, and the lens and the source, respectively. For con- 
venience, we will define a fiducial critical density for a 
source plane at z s = z ds = oo, and all models will be 
scaled to this standard. 

The convergence is related via a Poisson-like equation 
to a normalized potential: 

V 2 ^(0) = 2k(0). (3) 

Here and throughout this paper, all derivatives are in 
angular units in the lensing plane. A single light beam 
is deflected by: 

a = Vi/> • (4) 

The lens equation relates the source position (5 to the 
image position(s), 9, as: 



(3 = 6- Vip. 



(5) 



When the lensing potential does not vary appreciably 
across the source, the lens mapping can be linearized. 
The transformation between the source and the image is 
given by the Jacobian matrix 



1 - k - 7i -72 
-72 1 - K + 7l 



(6) 



From this, we see that distortions in shape are well 
described in terms of shear which is related to the lensing 
potentials through the relations: 



1 



7i = 2^' n _ ^' 22 ) 

72 = "0,12 



(7) 

(8) 



using Einstein convention for derivatives. 

The radial eigenvalue is given by, A + = 1 — k + [7I and 
the tangential eigenvalue is given by, A_ = 1 — K — I7 
. The matrix is singular where X± = 0. These points 
define the critical curves of the lens. 

Third order corrections to the lensing potential be- 
comes non-negligible when the lensing potential varies 
across the image. These obser yables, the gravi- 
tation al flexion, were derived in iGoldberg fc Bacon] 
(|2005h . They are more colloquiall y referred to as the 
"banananess" (Schneider & Er 20071) or bending of an im- 
age. It is based on the thir d angular derivatives of the 
potential (|Bacon et al.ll2006h given by, 



^"=(71,1 +72,2)+* (72,1 -71,2) (9) 

£/ = (7i,i -72,2) + i (72,1 + 72,2) (10) 

Each of the deflection (a), the shear, (71, 72), the con- 
vergence (re) and the flexion (T,G), are linear functions 
of the potential field. While the discussion in this pa- 
per primarily focuses on measurements with sources at 
fixed, infinite redshift, it should be noted that each of 
these terms scales as: 



and where 



k(z s ) = Z(z d ,z s )n(z s = 00) . (11) 



Z(zs) = ^ • (12) 



3. GRID-BASED CLUSTER LENSING 

Mass reconstructio n studies have been very s uc cessfu l 
on cluster scales (jBartelmar m & Schneider] (|2001f ); 
I Clowe fc Schneider! (12001): iHoekstra et alJ (|2002t i: 
iBroadhurst et all (|2005bD : IOkura et all (|2007h and refer- 
ences therein). Because these systems typically contain 
many lensed images, the shear signal can be extracted 
with high significance. In this section, we describe an 
important class of cluster inversion techniques which 
reproduce the convergence field on a grid. Our specific 
choices of grid-based techniques include those which 
have already been extended to include strong-lensing 
information with non-parametric models and thus 
provide a fertile basis for comparison. Further, there 
are many variants even within the sub-category of grid- 
based reconstruction techniques. We focus primarily on 
th eir commonali ti es, as exempli fied by those d is cusse d 
in iBradac et al.l (|2005al |bT) and ICacciato et al.l (|2006t ) . 
We focus on methods in which various scalar fields 
{i/j, re} are defined on a Cartesian grid, and minimized 
according to the criteria described below. In so doing , 
we note some in t eresti ng exceptions: iDiego et al.1 (2005) 
and IS aha et al.l (|2001l ), who describe an adaptive mesh 
technique for refining the field o n differen t reso lution 
scales and iMarshall et~ai1 (|2002h ; iMarshalll (|2006l) who 
use a variable smoothing scale for their weak lensing 
mass reconstruction. 

3.1. Weak Lensing on Grids 

The standard app roach to lensing arclet inversion 
(jLuppino et al.lfl999h has been to measure the elliptic- 
ity of observed images as an unbiased estimator of the 
reduced shear: 

(e)=.9^ T J -- (13) 

1 — K 

For relatively weak fields (k <C 1), this is very nearly a 
direct estimate of the shear, and can perform a direct 
finite inversion to estimate the density field. 

In recent years, there has been a flurry of work on 
optimal me thods for non-param et ric cluster mass recon- 
struc tions (|Bradac et al.l l2005al fbl: iNatarajan fc Springell 
2004). In general, these papers focus on estimating the 
potential {ip} or convergence {re} fields of a cluster by a 
X 2 minimization analysis. Both the shear and the con- 
vergence are linear functions of the potential field. Thus, 
if a model potential field, {iJj}, is defined on a grid, then 



3 



the shear at some grid-cell, i, may be expressed as a lin- 
ear combination of potential: 

Tii = G$1>j (14) 

with a similar expression for the convergence, k, and the 
imaginary component of the shear, 72. We refer to these 
below simply as ^({ip}), since we wish to remind the 
reader that the estimate of the shear is an explicit func- 
tion of the test potential field. Because these fields are 
combinations of second derivatives of the potential field, 
the G^ 1 ) matrix and the others are easy to compute using 
finite differencing, and are extremely local. A very good 
graphical repres entation of the finite d ifference operators 
can be found in lBradac et~all (|2005bh . 

In the weak field limit, the complex ellipticity of a 
lensed galaxy is a linear, albeit noisy, estimator of the 
complex shear field. The principle component of noise is 
the intrinsic ellipticities of galaxies which follow a Gaus- 
sian distribution with standard deviation of a e ~ 0.3 for 
each component. The large variance in intrinsic elliptici- 
ties necessitates averaging over many images so that the 
noise in a single grid cell is zero or apply an artificial 
smoothing scale to a more finely gridded mesh. For a 
weak-lensing only calculation, a x 2 minimization is per- 
formed on: 

( 7,(W) \ 2 

xh^E ^^P J (15) 

i * 

where the estimate of n is taken from the previous iter- 
ation of the potential field, and thus, the model rapidly 
converges to a maximum likelihood solution to the po- 
tential field. 

3.2. Strong Lensing 

A n umber of researchers, including Brad ac et al.l 
(|2005al lbl) have noted that a similar grid-based formalism 
may be used with strong lensing signals. Strong+weak 
(S+W) reconstructions use both shear fields and the po- 
sitions of multiply imaged sources can be used to ac- 
curately reconstruct both the cores and halos of clus- 
ters. While our current PBL implementation, described 
in the next section, does not currently incorporate strong 
lensing analysis, we introduce this component of grid- 
based lensing reconstructions to illustrate how directly a 
strong- lensing analysis could be incorporated into PBL. 

Strong lensing by clusters produces an especially ele- 
gant result because if, say, two images are observed at 
positions, 6 , and 8 , then it must be true that both 
images originated at the same (unknown) position in the 
sky. Thus, we have a simple relation: 

6 A - a(6 A ) = e B - a(9 B ) (16) 

The appeal of this relationship is that it is fundamentally 
linear and thus the angular separation between the two 
images (itself, a measurable quantity), can be directly 
related to the difference in the first derivatives of the 
potential at two different points in the field. 
As above, the local derivatives can be computed as: 

a xi = Alf-ipj 

with a similar expression for the y component of the dis- 
placement. The matrix elements of A are easy to com- 
pute as they are simply the 1st derivative in a simple 



grid-based 2nd order difference scheme. More generally 
we can express this as ai({ip}). Thus, an additional x 2 
term can be added: 

2 _ (( a A (m)- a B (m))-(o A -o B )y 

xs- 2. -f 

impairs 

(17) 

and minimized either independently, or simultaneously 
with the weak lensing component. 

3.3. Regularization 

Using a \ 2 minimization technique discussed in § 13.11 
it is possible to get a checkerboard pattern due to inde- 
pendent noise in the two components of ellipticity. This 
requires the addition of a regularization term to the x 2 
to suppress this noise. 

Scale refinement is also necessary in cases in which 
strong+weak lensing signals are combined. To make 
this argument concrete consider a toy isothermal sphere 
model of a cluster with a 1-d velocity dispersion of 600 
km/s. Each multiply imaged pair will be separated by 
twice the Einstein radius, about 20 arc-seconds in this 
case. This represents the minimum necessary resolution 
in the reconstruction to say anything about strong lens- 
ing. 

On the other hand, even very efficient space-based 
weak lensing analysis of clusters seldom yield more than 
approximately 100 images/square arc- minute. Using a 
simple Poisson noise estimate, we may achieve uncer- 
tainty of er 7 = 0.06 only with images binned on scales 
larger than 30 arc-seconds on a side. Smaller binnings 
will naturally yield larger uncertainties. Simple grid 
based method cannot capture both the weak-lensing sig- 
nal to high accuracy as well as resolve the strong lensing 
regime. In order to deal with this issue, different inves- 
tigators have used different regularization techniques. 

One method is to use a series of finer and finer grid- 
dings, and at each successive level of refinement the 
convergence field from the previous l evel is matched as 
closely as possible. The lBradac et al.l (|2005b[) S+W tech- 
nique uses this method, with the weighting parameter 
selected to provide a \ 2 P er degree of freedom equal to 
1, such that: 

^ = ^EK (n) -^ 1} ) 2 - ( 18 ) 

i 

Where nf 1 ^ represents the estimated convergence on 

the previous, coarser, gridding, and where nf** 1 = 0. 
We use this form explicitly in |J5] where we test the 
PBL method and contrast it to grid-based reconstruc- 
tion methods. 

3.4. Some Questions 

Grid-based reconstructions have produced some excel- 
lent measurements, however, there remain a number of 
complications. First, grid-based techniques are really op- 
timized to measure a single scale, the grid-spacing. How- 
ever, as we discuss above, in many interesting systems, 
both the structure and information are hierarchical. An 
optimal technique should provide higher resolution in re- 
gions of greater information content. 
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Moreover, the smoothing and weighting of the strong 
lcnsing, weak lensing, and regularization are created in 
an ad hoc basis. The ideal smoothing scale should be 
variable, and such that the signal/noise ratio of the re- 
constructed field is similar in every smoothed cell. 

Third, the information from the image ellipticity can 
only be inverted outside the critical cur ves of the lenses. 
Insid e the (tangenti a l) cri ti cal curve (ISchneider et al.l 
1991 IPetters et al l l200lt ISchneider fc Weisd 1199a 
Hoekstra et al. 2004) there is an abrupt switch in par- 
ity of the induced ellipticity of an image. More plainly, 
in the regime I7I > |1 — k\, the ellipticity is related to the 
shear via: 



■rang 



g* 



(19) 



As discussed in £14.41 this produces a discontinuity in 
the ellipticity as a function of n and 7. No simple linear 
minimization scheme, even an iterative one, will converge 
to the "strong lens" solution if one starts with a "weak 
lens" initial guess for the local potential field. 

4. PARTICLE BASED LENSING - PBL 

In this section, we introduce a new technique called 
Particle Based Lensing (PBL) which has the ability to 
combine the disparate lensing scales in a coherent way 
without requiring a regularization scheme. Several of 
the concerns discussed in the previous sections have to 
do with the method of discretizing the data for the re- 
construction of the lens potential. In order to address 
this, we turn to a technique which is widely used in an- 
other area of astrophysics in which information must be 
analyzed on a wide range of physical scales - numerical 
N-body simulati ons. Smoothed Pa rticle Hydrodynamics 
(SPH; see, e.g. iMonaghanl (|2005l h for a recent review) 
is used in the modeling of a wide range of physical sys- 
tems including planets (Woolfso n1l2007f). star form ation 
(|Springel fc Hernquistj l2003t iNagamine et all l2004h and 
galaxy formation ( Kaufmann et al.ll2007l ). The math- 
ematical details of PBL can be complicated, hence we 
have made our codes for the method public 1 through our 
website. 

Before getting into the details, however, it is impor- 
tant to emphasize what PBL is and is not. PBL is a 
new way of discretizing and describing a reconstructed 
field. Moreover, it includes a metric for comparing a re- 
constructed model to the observed data. Everything we 
describe below is aimed at demonstrating why this model 
and metric are ideal for lensing systems with uneven 
information content. While the current code, and the 
worked examples are based on weak-lensing data only, 
PBL is based on the idea that other probes of the poten- 
tial field: strong-lcnsing positions, flux ratios, and flex- 
ion, can be added to the metric with little complication. 

PBL is not, however, a minimization scheme. That is, 
much like grid-based reconstruction methods, PBL fun- 
damentally consists of a list of dimensionless potentials 
and a metric to describe the goodness of fit. It does not 
describe how that minimization criterion is to be met, 
however. In our model section, we describe a number 
of approaches to efficient model convergence. The ma- 
jor argument in favor of PBL, however, is not that \ 2 

1 http : //www . physics . drexel . edu/~deb/PBL . htm 



minimizes efficiently, but rather that a low \ 2 in PBL 
actually corresponds to a model which closely matches 
the true underlying system. 

4.1. A Particle Description of the fields 

The fundamental description of the PBL field lies in 
the a list of potentials, {ip}, one each at the positions 
of each observed lensed image. In order to make the 
field as continuous as possible, we may expand the local 
potential field around the position of any lensed image, 
(tpn, in this case) to arbitrary order: 

lp(0) = Vn + Ojlpnj + ^ W A<. + ... (20) 

where is the relative offset of the test-point from galaxy 
n. 

As with grid-based lensing, the local derivatives are 
composed of a linear combination of the potentials at 
each grid point. That is: 



■111 



(21) 
(22) 



and so on for arbitrarily higher derivatives. In reality, we 
typically extend the D^ u > matrices up to 3rd order, where 
v corresponds to 2 matrices for 1st derivatives (displace- 
ment field), 3 for second derivatives (shear and conver- 
gence), 4 for 3rd derivatives (flexion). Here we use Ein- 
stein summation convention, the sum over m runs from 

1 to Ng. 

In terms of the D^ u > matrices, equation (j2"0|) may be 
rewritten as: 



(23) 



where we are explicitly estimating the potential at the m- 
th galaxy from the local derivatives defined at the n-th. 
We also compactify equation (|20[) by defining: 



x {1) -e - 

y(2) _n 

X (3) =- (6 



ray 



(24) 
(25) 

(26) 



and so on. 

In order to estimate the derivatives of the potential 
field near each galaxy, we need to first compute the 
matrices. Since this problem is under-constrained, we 
solve for these matrices via a \ 2 minimization: 



^m-^n-y^D^X^A W nm (27) 



where w nm is a window function, guaranteeing that only 
neighboring galaxies effect the potentials of one another. 
We use a window function of the form: 



\0 n — Q-n 



(28) 



where h n is inversely proportional to the signal-to-noise 
at the n-th image positions. The smoothing scale can also 
be chosen to be of the form h nm , i.e symmetric between 
between the points n and m. 
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The signal-to-noise is a function of the local density of 
background images and type of constraint (e.g. elliptic- 
ity positions of multiple images, etc). A similar approach 
of using signal to noise dependent smoothin g scale has 
been u sed in image analysis of X-ray data lEbeling et all 
(|20Mk In regions where there is a high density of infor- 
mation, the smoothing scale h n may be set much lower 
than in regions of low information density. 

This function must be minimized for every matrix ele- 
ment such that: 



dx 



yv nm w nm I Ym Yn / J XJ ny yy nm'rP 



3D 

(29) 

But since equation (j29|) is under-constrained, we may 
also say: 

d 2 x 2 

(30) 



t) 







yielding: 



(31) 



for all n,m and v. This can be solved with a simple ma- 
trix inversion, yielding the desired elements for D^"\ Of 
course, since the elements are a function only of the posi- 
tions and weightings of the galaxy images, these elements 
need only be computed once. The method potentially in- 
corporates higher-order derivatives of the potential, thus, 
combination of strong, weak and flexion information be- 
comes a relatively straightforward minimization problem. 

4.2. PBL vs. Regularization 

One of the major advantages of PBL is that we no 
longer need to introduce an explicit regularization in or- 
der to resolve multi-scale structure in a reconstruction. 
The various regularization schemes discussed in § 13.31 are 
not motivated from the associated observations, but are 
rather derived from assumptions about the mass profile 
of a cluster motivated by theory and simulations. 

However, one of the motivations behind using gravi- 
tational lensing is to be able to measure the projected 
mass without making any assumptions about the physi- 
cal state of the system. The advantage of using PBL is 
that we do not need to make any assumptions that go into 
choosing the regularization term. The smoothing scale 
of a "pebble" is controlled by h n which is determined by 
the local signal to noise. This means that the position 
representing weak lensing measurement will have a low 
signal to noise and correspondingly a high h n . This is 
similar to the typical weak lensing measurement which is 
done by averaging over a bin size larger than ~ 30". In 
case of strong lensing we know the positions of the multi- 
ple images for certain, implying high signal to noise and 
correspondingly low h n . This can be a few arc-seconds 
which is the scale at which the strong lensing structure 
can be resolved from multiple images. Thus scales of a 
few arc-seconds can be combined with scales greater than 
~ 30" without making any assumptions about the mass 
profile, rather by taking input from the data. 

4.3. Estimation of the Potential Field 



As with grid-based lensing analysis, in PBL, we use a 
X 2 minimization to estimate a maximum-likelihood po- 
tential field. In this case, however, we sample the poten- 
tial at every point, and use the local derivatives of the 
potentials as defined in equation (|22p to minimize: 



i=2\ n=N„ 



X' 



E 



i 



7i t) (M) 



r(0 



1 



(32) 



where i ranges from 1,2, and indicates the real or imagi- 
nary component of the shear, reduced shear, or ellipticity. 
=W8: shall henceforth refer to the first term in the paren- 
theses as gn \{ip}), the estimate of the reduced shear of 
a model, and the weighting term outside the parentheses 
as w n , yielding: 



X 



(33) 



which is the form we will refer to from now on. 
This is a weak lensing only expression. Replacing 

9n\{ip}) with l/gn({ip}) gives the strong lensing coun- 
terpart of Eq. [3S] In the next section we discuss how we 
include this strong lensing version of the equation. 

4.4. Interpolated Ellipticities 

Linear inversion techniques require that the function 
to be minimized is smoothly varying over the domain of 
interest. The ellipticities are given by two functions in 
the weak and strong lensing regimes by Eqs. (|13ll9p . The 
boundary of the two regimes define the critical curves 
where \g\ — 1 making ellipticities continuous but not 
diffcrentiable. 

The transition between the two regimes can be fa- 
cilitated if the sources are distributed in redshift, but 
minimization functions will be much easier if we allow 
a smoothing of the discontinuities. This is a two step 
process, first we need to write Eq.. (|13|19|) in terms of a 
step function, 



[!-«(<?)]<?• 



n(g) 
g* 



(34) 



where the function TL(g) is a step function at g = 1. We 
may replace the step function by an approximate smooth 
function. We define: 



(35) 



Here 770 is the free parameter that controls the accuracy 
of the step function. A higher value of ijq makes the step 
function more accurate. The step function is approxi- 
mated as (Fig. [TJ, 



H{u) = 



1. 



1 



-2 a 



(36) 



This approximation replaces the ellipticities only in the 
neighborhood of the critical curves (discontinuity) by a 
continuously diffcrentiable function. The problem can 
now be solved by standard minimization techniques. The 
interpolated ellipticity function is shown by a dotted line 
in the second panel of Fig. [TJ showing the derivative dis- 
continuity explicitly. 
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Fig. 1. — In the upper panel, we plot the interpolated Heavi- 
side step function. It is clear from the plot that the function is 
only approximated by a smooth function near g=l, for all other 
g it behaves like an ordinary step function. Also higher value of 
the parameter r)o increases the accuracy. In the lower panel, we 
plot the resulting ellipticity as a function of reduced shear for the 
combination, \^\ = n. 

4.5. x 2 Minimization 

When we first introduced PBL above, we remarked 
that it was primarily a way of describing a lens recon- 
struction in such a way that a small \ 2 would necessarily 
correspond to a good representation of the underlying 
field. In practical terms, though, for a reconstruction 
code to be useful, we need to describe a means of mini- 
mizing (or nearly minimizing) the x 2 . Below, we describe 
our pipeline for fast convergence of a maximum likelihood 
solution. 

While PBL is a non-parametric reconstruction scheme, 
it has the useful property that we may start a minimiza- 
tion with any assumed model we like. However, no extra 
weight is given to our a priori assumptions. At the end 
of a minimization we may simply use the standard tech- 
niques to estimate the likelihood of a particular value of 
X 2 - 

That said, even with the caveat above regarding 



smoothing of critical curves, it is very difficult to 
smoothly vary a solution such that stron gly lensed re- 
ions are produced. As pointed out by iBradac et al.l 
2006) a x 2 minimization process does not ensure reach- 
ing a global minimum. 

To that end, our initial configuration of {tp} is gener- 
ated by laying down a small number of Singular Isother- 
mal Spheres (SIS's). Since there are a low number of 
parameters (3 for each model sphere) , a global minimum 
may be reached through a combination of trial and er- 
ror, simulated annealing, or even (for small numbers of 
spheres), finite sampling. Indeed, one may even use an 
interpolation of a reconstruction recommended by a grid- 
based solution. For systems with strong lenses, one may 
apply the recons tructed field generated by "LensPerfect" 
(|Coe et al.ll2008| ). for example as a starting point. 

We hasten to remind the reader that while this tech- 
nique will produce the optimum parametric fit, it will 
not, in general, produce the overall best fit. As a result, 
further iteration is required. 

We have found that by starting with an initial model 
with well-identified strong- lensing regions, convergence 
to x 2 1 DOF ~ 1 may be achieved relatively quickly, even 
if the strong lensing regions are only approximate. For 
the current implementation of our code, we use New- 
ton's method to reach a local minimum. We have found 
satisfactory, fast, convergence for several thousand back- 
ground sources. 

5. TEST APPLICATIONS 

In this section, we apply PBL to three systems as 
a proof of concept. In the first, we model a Softened 
Isothermal Sphere, and examine the relative abilities of 
PBL and grid-based inversion to reconstruct the a rel- 
atively peaked core. In the second, we model a super- 
position of two softened isothermal spheres at a given 
separation as a simple model of a system with sub- 
structure. Fi nally, we reconstruct the "Bullet Cluster" 
(1E0 657-561 (iMarkevitch et all l2002t 12004 IClowe et all 
jOflSlBradal et alJl2006b iClowe et al.ll200l . an observed 
multi-peak system of considerable interest. We show that 
using weak lensing alone, we are able to reconstruct both 
Dark Matter peaks. 

5.1. Simulation: Softened Isothermal Sphere 
5.1.1. Model 

We begin by generating a softened isothermal sphere 
with a potential: 



and convergence: 



1> = 0eVO^+¥ c , 



(6 2 + 29 2 c ) 



\3/2" 



(37) 



(38) 



where 6e is the Einstein deflection angle given by 

The data is simulated on a unit square field of view. 
For simplicity we have assumed all sources to be at 
z = oo, with E = 0.2 and 9 C = 0.08. We lens 607 back- 
ground galaxies, and apply an intrinsic ellipticity (noise) 
with o~ Ea = 0.1 in each of the principle directions. For 
all further calculations we use a ACDM cosmology with 
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f2 m = 0.27 and £l\ = 0.73. This configuration represents 
a galaxy cluster at a redshift of zi ens = 0.4 with a veloc- 
ity dispersion of a v — 850 km/sec. The field of view is 
105" and 0.5 Mpc. 

5.1.2. PBL and Grid-Based Reconstructions 

For the single peak and the double peak simulation(see 
below), we perform both a grid-based reconstruction as 
well as PBL . We us e the regularization suggested by 
iBradac et al.1 (|2005bf ). and described in detail in § 13.31 
for the grid based method. In case of the single peak 
the reconstruction is initially performed on a coarse grid 
[n x = 6 gridcells), and is refined up to n x = 24, using 
the k estimated at each previous step as the prior. For 
the double peak system we start with n x = 10 and refine 
up to nx=40. For both systems the final reconstruction 
contains less than one particle per grid cell. 

For the PBL reconstruction, we use a smoothing scale 
of the form: 

K = 1 ~T7 (39) 

where p n is the local number density of points , c is 
a constant, and £ is a tunable parameter to maximize 
signal-to-noise. For our simulation £ = 1 is an optimal 
choice and for the observational case we have used £ = 
0.5, which is the obvious choice for equalizing signal to 
noise for every smoothing length. We select c such that 
the integrated signal-to-noise is greater than unity. The 
PBL reconstructions are gridded to the same resolution 
as the grid-based reconstruction to aid visualization. 

For both reconstructions, we begin our iterations with 
a best-fit SIS. We do not, however, use this in the regu- 
larization for the grid-based reconstruction. 

5.1.3. Results 

Before discussing the results, we note a potential com- 
plication. Gravitational lensing mass measurement suf- 
fer from mass sheet degeneracy when the sources are not 
distributed in redshift. This implies that k can be deter- 
mined up to a degeneracy Xk + (1 — A). This transforms 
to a degeneracy in the potential of the form, 

4>(8)^^(e) = ^{\-\)e 2 + \^{e) (40) 

For the simulated data we have computed the best value 
of A in each case and transformed our reconstruction with 
that value of A for both the grid based method and PBL. 

In Table [1] we compare the x 2 for the best fits of both 
the grid-based reconstruction along with PBL for a vari- 
ety of smoothing normalization parameters, c. The aim 
of this table is to quantify the deviation of the recon- 
structed k from the true k. In each case, the ostensible 
X 2 /DOF is of order unity. However, one needs to be 
careful with simply asserting that the lower x 2 produces 
the best result, since the regularization in grid-based re- 
construction adds a penalty function, and the smoothing 
scale in PBL lowers the effective degrees of freedom. 

So while both models produce small values of x 2 > the 
real question is whether these good fits correspond to 
a an accurate reconstruction of the underlying density 
field. In Table [l] we do several comparisons which re- 
late the reconstructed k at each galaxy (or grid-point) 
with the true k modeled by the simulation for both the 



single peak and the double peak. The comparisons are 
done with a range of values for both rj (the regularization 
weight in grid based method) and c (the proportionality 
constant in PBL). 

The first column in Table 1 describes the method used, 
i.e either PBL or the grid based method. The second 
column describes the difference between reconstructed k 
and the true k for every galaxy position. In order to ex- 
tract this information from the gridded reconstruction 
we have used the nearest grid point method which sim- 
ply means the k at each galaxy position is assigned the 
value at the corresponding grid cell. The third column 
describes the deviation of the reconstructed k from the 
true k at every grid cell weighted by the number of im- 
age galaxies in that grid cell. The 4th column describes 
the difference between reconstructed and true k weighted 
uniformly over the grids. In each of the 3 compar- 
isons, PBL reproduces the original reconstruction with 
the highest fidelity. The 5th column gives the x 2 /DOF, 
the 6th gives the regularization parameter rj for grid 
based method and the 7th column gives the smoothing 
normalization parameter for PBL. 

In Fig. [21 we show the radial reconstruction of the 
softened isothermal sphere using the two different tech- 
niques. The bulk of the penalty associated with the 
grid-based reconstruction relative to PBL occurs near the 
core. By construction, PBL is designed to perform well 
in this regime. 

5.2. Simulation: A Double Peaked Cluster 

5.2.1. Model 

While PBL has been shown to perform well model- 
ing a single Softened Isothermal Sphere in the previous 
section, the other major goal of this method is to re- 
construct small-scale substructure in a system. To that 
end, we model a doubly-peaked system with 814 lensed 
background galaxies. As before, they are placed on a 
unity grid, and are modeled as 2 Softened Isothermal 
Spheres, with: x\ = 0.65, y\ — 0.35, x 2 = 0.35, y% = 0.65, 
E i = 0.2, 9 cl = 0.1, 6 E2 = 0.2, and c2 = 0.1. The simu- 
lated noise, and reconstruction technique for the double 
peaked system are identical to the single peak system. 
This is system of two sub-clusters at a redshift of zi ens = 
0.4 having a velocity dispersion of o~ v = 850 km/ sec sep- 
arated by 226 kpc. The field of view is 105". 

5.2.2. Results 

As with a single sphere, both PBL and grid-based re- 
constructions produce x 2 /DOF ~ 1, as illustrated in Ta- 
ble [1] However, as with the single sphere reconstruction, 
PBL produces smaller errors with regards to the under- 
lying model than does the grid-based reconstruction. 

In Fig. [31 we show a grey-scale plot of the residuals be- 
tween the underlying model and each of the reconstruc- 
tions. Unsurprisingly, both models have the greatest dif- 
ficulty reproducing highly peaked cores, though PBL is 
more responsive to high local gradients in k. We describe 
the general quality of the fit in Table [T] 

5.3. Observation: The Bullet Cluster 

5.3.1. Observations 

Finally we perform a mass reconstruction of the bul- 
let cluster (1E0657-56). This galaxy cluster is a rare 
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TABLE 1 

Comparison between PBL and grid based method." 



Method ^ — 


K model,i) 2 


N 2 ., „ N 2 
Ej 9 ("i-KmodeM) ™« E,» r 




X 2 /DOF 










N 2 ., 






V 


c 






v^i grid 


N grid 








Single Peak 














PBL 


0.0200 


0.0136 


0.0147 


1.03 




0.5 


PBL 


0.0181 


0.0128 


0.0119 


0.94 




0.7 


PBL 


0.0219 


0.0139 


0.0131 


0.95 




1.0 


PBL 


0.0235 


0.0140 


0.0133 


0.94 




1.3 


PBL 


0.0227 


0.0120 


0.0121 


0.98 




1.5 


GRID 


0.0311 


0.0283 


0.0237 


0.6 


10 




GRID 


0.0309 


0.0280 


0.0223 


0.79 


30 




GRID 


0.0311 


0.0280 


0.0224 


0.94 


60 




Double Peak 


PBL 


0.0250 


0.0174 


0.0167 


0.82 




1.1 


PBL 


0.0231 


0.0168 


0.0160 


0.80 




1.38 


PBL 


0.0277 


0.0193 


0.0180 


0.82 




1.7 


PBL 


0.0320 


0.0219 


0.0208 


0.87 




2.0 


GRID 


0.0570 


0.0711 


0.0630 


0.92 


20 




GRID 


0.0367 


0.049 


0.039 


0.7 


10 




GRID 


0.0359 


0.0482 


0.0454 


0.83 


60 





a The 2nd, 3rd and the 4th columns represent the deviation of the reconstructed re from the true re weighted 
uniformly by galaxy, by local density within gridcells, and uniformly by gridcells. Here r\ is the weight given 
to the regularization for grid based method and c is the smoothing normalization parameter for PBL. 



CD 



0.0 0.2 0.4 0.6 0. 

r 



0.2 0.4 
r 



0.6 



Fig. 2. — A radial plot of the reconstructed convergence (re) of a simulated Softened Isothermal Sphere. The circles r epre sent binned 
reconstructed re and the error bars represent the scatter in each bin. The dots represent the true value of re given by Eq. 1381 The x-axis 
represents the radial distance from the center of the softened isothermal sphere. The radial distance is scaled and hence unit-less. Upper 
Panel: Using PBL. Lower panel: Using grid based method. The error bars in the radial plot using PBL is higher. This is because the 
errors introduced in PBL are dependent on the local signal to noise which are not spherically symmetric. In the grid based method the 
errors are averaged uniformly on the length scale of a single grid which makes the radial scatter very low. 



supersonic merger in the plane of the sky. Its distinc- 
tive structure and orientation makes it an ideal cluster 
for observing dark matter using gravitational lensing. It 
consists of two sub-clusters separated by 0.72 Mpc, which 
have just undergone a merger and are moving away from 
each other. The western sub-cluster is less massive and 
the eastern main cluster is more massive. The line-of- 
sight velocity difference suggest that their cores passed 
each other 100 Myr ago. The collisionless dark mat- 
ter in each of the sub-clusters have crossed each other 
but the fluid-like intracluster plasma is in the process of 
electromagnetic and thermal interaction producing high 



X-ray luminosity far removed from le nsing mass peaks 
(|Clowe et al.ll2006t iBradac et al.ll2006f ). 

For the Bullet Cluster, we perform a PBL reconstruc- 
tion only, since it has been well -s tudied with g rid-based 
methods (using ISchneiderl (| 1 995ft i lKalserl (|1995| i) and the 
K-contours are publicly available. We use publicly avail- 
able weak lensing data from the Bullet Cluster Project 
Page 2 . The catalog was constructed using data from 
three different instruments: the ESO/MPG Wide Field 
Imager, IMACS on Magellan, and two pointings of ACS 

2 http: //flamingos . astro .ufl . edu/le0657/public .html 
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40 



Fig. 3. — The plot of the difference between reconstructed 
convergence^ and true re for the double peak SIS system. Left 
Panel: Using PBL . Right panel: Using grid based method de- 
scribed in § [3] Both maps are gridded for easy visualization. Also 
there are empty grid cells with no image galaxies. The value for 
those grid cells in the above difference map is set to zero for both 
reconstructions. As we can see the error in the cores of the peaks 
is much lower using PBL mass reconstruction 

on HST. The shapes of the galaxies were measured in- 
dependently on each of the image sets averaging for the 
common galaxies. The weighting for each galaxy is based 
on its significance of dete ction in every im age set and 
normalized appropriately (|Clowe et al.ll2006f h 

The catalogs were combined using weighted average re- 
duced shear measurements and the weights of individual 
galaxies were increased when they occurred in several 
catalogs. This weighting is listed in the shear catalog. 
We include this weighting in our reconstructions as well 
and choose only those images with a weighting greater 
than 1. As we have already illustrated in the simulations 
PBL is most effective when the information density is 



variable, i.e close to the core of the clusters. In case 
of the bullet cluster we zoom into a region bounded by 
104.53° to 104.69° in right ascension and 55.92° to 55.97° 
in declination. Following this cut, our sample includes 
1259 weak lensing background galaxies. In order to do 
the mass reconstruction we use the average redshift of 
this sample, z = 0.91. 

5.3.2. Reconstruction 

The Bullet Cluster was made famous by the direct de- 
tection of dark matter by IClowe et al.1 {2006). Indeed, 
since one of the major findings of this group is that the 
dark matter appears offset from X-ray emissions, we do 
not include any prior model when reconstructing the sys- 
tem, but are able to achieve fast convergence with two 
clearly visible peaks. This reconstruction guides us in 
choosing an initial condition for subsequent \ 2 minimiza- 
tion. 

We have calculated the integrated mass within 150 
kpc of each peak. The main peak has a mass of 1.57 x 
10 14 M(D and the su b-cluster has a mass of 0.9 x 1O 14 M . 
IClowe et all (|2004D report a value of (1.02 ± 0.16) x 
10 14 M Q for the main peak and (0.66±0.19) x 10 14 M Q for 
the sub-cluster within 150 kpc of the each peak. In each 
case, our estimate exceeds that of Clowe et al. by approx- 
imately 3.4 er. However, more a more recent S+W recon- 
struction by the same group (jBradac et al.l 120061 ) yields 
masses of (2.8±0.2) x 10 14 M Q around the main peak and 
(2.3 ± 0.2) x 1O 14 M around the sub-cluster within 250 
kpc of each peak. Inclusion of strong lensing information 
makes reconstruction of the cores more accurate and also 
leads to a higher estimates of the mass. Even correcting 
for the greater area, this suggests Clowe's initial mass 
estimate may have been low. 

Our mass estimates usi ng PBL is hig h er tha n the weak 
lensing reconstruction of IClowe et alj (|2004f ). and thus 
more in line with the S+W results. This is a result of a 
difference in method. For example, we start from an ini- 
tial condition and iterate to the correct solution whereas 
IClowe et al.l (|2004T ) have fitted a radially averaged shear 
profile to the NFW or King profile. As already seen in 
the simulations using an initial condition recovers values 
of k close to the core with greater accuracy. This implies 
that while most weak lensing k maps report K-contours 
less than 1 using initial condition and PBL we are able 
to get k greater than 1. This implies that the mass we 
measure will also be greater than the typical weak lens- 
ing mass measurement. Also t o measure the mass of the 
sub-cluster lClowe et al.l (|2004f ) have removed the mass of 
the main cluster to avoid over-estimation of the mass, we 
have not considered this effect in our reconstruction. 

In Fig. [4j we show our PBL reconstruction of the bul- 
let cluster. Note that, despite using weak lensing signals 
only, we are able to identify both density peaks and using 
initial conditions we are able to get k > 1 for the main 
peak. We also do a comparison of the publicly available 
K-contour with the K-contours reconstructed using PBL. 
The location of the main peak coincides for both recon- 
struction. The sub-cluster contours for PBL are slightly 
removed from the publicly available K-contours. 

Error analysis for PBL will be discussed in detail 
in future papers. In particular the noise covariancc 
matrix, ((n — (k))(k — (k)) t ), will give us important 
insights into the errors caused by the reconstruction 
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Fig. 4. — A weak-lensing only reconstruction of the bullet cluster 
using PBL described in § [3] Note that both substructure peaks are 
clearly identified. Upper Panel: This the K-map using PBL. The 
cross denotes the centroid of the multiply imaged positions. Lower 
Panel: This a comparison of the n contour derived using PBL(solid) 
and the publicly available contour plot of ft(dashed). 

method. A bootstrap method can also be used to de- 
termine error bars on mass measurements from observa- 
tions. In case of simulations several Monte Carlo realiza- 
tions of the noise can be used to study the errors. 

6. DISCUSSION AND FUTURE PROSPECTS 
6.1. Additional Signals 

Thus far, we have developed the formalism for PBL, 
and done worked examples demonstrating how it may 
be applied to weak-lensing reconstruction. It is designed 
to model structure hierarchically, in part because of the 
great success of Strong+Weak lensing analysis. 

Several groups have already shown how multiple im- 
age positions may be added to the information yielded 
by lens ellipticities to produce very high quality mass- 
maps of clusters. It was our desire to maximally exploit 
the different information scales of the strong and weak 
lensing signals which motivated the development of PBL 
in the first place. 

However, there is yet more information besides image 
differences potentially available which may be utilized 
in a reconstruction. Consider that in addition to the 
two constraints generated by the positional difference be- 
tween two images, we also can measure a flux ratio, and 2 
ellipticity differences. Thus, in principle, we have 5 mea- 
surable, model parameters per strong lensing pair rather 
than 2, and in an idealized case, this improves potential 
resolution of a system in the strong lensing regime by 
V(5/2) -1.6. 

As a way of guiding the future development of PBL, 
we discuss possible future avenues of investigation below. 



6.1.1. Flux 

Apart from the centroid position, the Petrosian flux of 
an image is the most straightforward to measure. The 
relationship between magnification lens is simply the in- 
verse of the determinant of the projection matrix: 

M = (1-k) 2 -| 7 | 2 (41) 

Unlike the displacement vectors (a), which are simple 
linear operators of the potential field (the gradient), or 
the weak-lensing shear field which is nearly so (since in 
the limit of k << 1, the image ellipticity is an unbiased 
estimator of the shear field) , the flux is a highly nonlin- 
ear function of the shear and convergence fields. This 
accounts, in part, for the reason that it has not been 
used previously i n clust er reconstructions. Here we note 
that ISaha et all (|2007l ) show that the image positions 
itself constrain the fluxes for a source with three non- 
collinear components. This is a special case, for cluster 
lensing three component sourc es for strong lens i ng may 
not always be available. Also, iNatarajan et al.l (|2007al ) 
use magnification information in their parametric mass 
modelling of clusters. 

The other major consideration is that magnification is 
not a smoothly varying function of the potential fields. 
It is well-known that on th e critical curves, magni fication 
goes to infinity (see, e.g. ISchneider et al.l (|1992| ) for an 
extensive discussion), but this is a set of measure zero, 
so in and of itself produces no problem. The issue is that 
the parity of the image reverses as an image that crosses 
the critical curve. 

Negative magnification means nothing more than re- 
versal of image parity, and thus cannot generally be easily 
detected. Thus, we are much more interested in comput- 
ing terms which scale like . Indeed, since we cannot 
measure the magnification directly, but only the flux, we 
propose that the combination: 



/I 



fi + fl 



(42) 



is directly measurable, and has no poles. 

Even so, a lensing model predicts a parity for a partic- 
ular image, and as with ellipticity, minimization, there 
is a discontinuity in the derivatives. In Fig. [5] we show 
the magnification (including sign) as a function of con- 
vergence and shear. 

6.1.2. Ellipticity Differences 

Likewise, while most measurements of the shear are 
based on an assumption that any given image is ran- 
domly oriented, two images of the same source are not. 
The difference in their measured ellipticity can be wholly 
modeled by the relative lensing fields at their respective 
locations. If both images were in the weak regime, we 
would be able to use the simple estimator 



£A - £B = ^ 1A - IB 



(43) 

where all terms in the equation are complex, and thus 
provide two constraints with high signal to noise per im- 
age pair. 

In general, however, a more likely configuration is that 
one image may be in the strong regime, and one in the 
weak. If we can determine from the configuration of 
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Fig. 5. — The magnification as a function of shear and conver- 
gence. The lower panel is a simple slice through the upper, with 
the choice 7 = k. Neither the magnification nor its derivatives are 
a continuous function. Moreover, flux ratios are only measurable 
for systems with at least two images (obviously). One or more 
of the images will necessarily have negative parity. Thus, a solu- 
tion to the potential field which is found using standard relaxation 
methods will not normally converge to a negative parity estimate 
for any magnification. 

lenses which is which, we might imagine a better esti- 
mator as: 

£A - £b = — - 9B (44) 
9a 

with the only associated noise corresponding to photon 
noise rather than random variance in the intrinsic ellip- 
ticity of the images. 

6.1.3. Flexion 

Thus far, the analysis of clusters in the weak or semi- 
weak re gime has primarily relied on shear. However , 
recently. lOkura et alj ([2007), and iLeonard et alj (|2007D 
have worked on reconstructing A1689 using flexion. In 
particular , the Okura group used a Fourier inversion sug- 
gested by iSchneider fc Erl (|2007f ). However, the advan- 



tage of our proposed PBL is that flexion (and, in princi- 
ple, any higher-order derivative of the potential) may be 
explicitly included as additional constraints in the clus- 
ter reconstruction. Unlike Fourier techniques, which rely 
on binning of the data, the PBL method will allow us to 
exploit the natural small-scale signal probed by flexion. 

6.2. Summary 

We have developed PBL, a new particle based tech- 
nique of mass reconstruction of clusters. The distinguish- 
ing feature of PBL is its ability to adjust its smoothing 
scale depending on the local signal to noise or the type of 
constraint and thus not require any regularization. PBL 
has the scope of calculating derivatives up to any or- 
der. Hence, lensing constraints that are a function of the 
derivatives of the potential can be easily included in the 
reconstruction. In this paper we have successfully ap- 
plied PBL to do weak lensing only mass reconstruction 
for a single peak and a double peak system. We have 
made the codes for PBL publicly available for applica- 
tion weak lensing measurements through our website(see 
§ |3| . The codes have been tested on the data sets and 
simulations described in the paper. A larger data sample 
will require modification of the current version of codes. 

As already explained PBL is a method of discretizing 
data and not a minimization method. A \ 2 minimiza- 
tion does not necessarily ensure reaching a global min- 
imum. In many cases the global minimum is guarded 
by steep walls surrounded by shallow valleys. Without 
any prior knowledge of the mass distribution it is very 
easy to get trapped in a shallow valley and not reach the 
global minimum. We have started with an initial condi- 
tion and interpolated the ellipticity function to aid us in 
this regard. 

In future work we will be including the additional con- 
straints, like the flux ratios, ellipticity differences and 
flexion along with measured ellipticitics and strong lens- 
ing positions. We will also be exploring different min- 
imization schemes to facilitate convergence to a global 
minimum. 
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